PVA of the NBI: Demography - GLMM of the number of fledglings per potential mother and per nest
- Hypothesis: We hypothesize that at the present time the northern bald ibis population can survive without further management and release. We predict that the observed demographic rates will ensure population growth and do not differ between the breeding colonies.
- Study area: Austria, Germany and Italy
- Data: Data of the number of fledglings per potential mother and per nest
We create GLMMs to investigate whether the raising type or colony has an influence on a) the number of fledglings per potential mother (incl. breeding probability) and b) the number of fledglings per nest (only actual breeding events and incl. temporarily added females as mothers).
Setup
## for non-CRAN packages please keep install instruction
## but commented so it is not run each time, e.g.
# devtools::install_github("EcoDynIZW/template")
## libraries used in this script
## please add ALL LIBRARIES NEEDED HERE
## please remove libraries from the list that are not needed anymore
## at a later stage
library("DHARMa")## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car
## This is DHARMa 0.3.3.0. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa') Note: Syntax of plotResiduals has changed in 0.3.0, see ?plotResiduals for details
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
Data
f_fledg_mother <- readRDS(file = here("output", "data-proc", "f_fledg_mother.rds")) # from script "05a_Reproduction_analysis" # chicks_year
f_fledg_mother <- droplevels(f_fledg_mother)
nest_numbers_fm <- readRDS(file = here("output", "data-proc", "nest_numbers_fm.rds")) # from script "04_Survival_analysis" # NBI_breeding_all
nest_numbers_fm <- droplevels(nest_numbers_fm)Female fledglings per potential mother per year
Based on the values used for RR Baseline
## [1] "parent_f" "hatch_year" "nr_f_fledg" "raising_type" "breeding_area"
GLMM of the number of fledglings per potential mother per raising type
Data exploration
# Plot the number of fledglings
ggplot(data = f_fledg_mother, mapping = aes(x = nr_f_fledg)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Plot of the number of fledglings per raising type per year
f_fledg_mother %>%
group_by(raising_type, hatch_year) %>%
mutate(sum.fledg = sum(nr_f_fledg)) %>%
ggplot(mapping = aes(x = hatch_year, y = sum.fledg)) +
geom_line(aes(group = raising_type, color = raising_type))##
## foster parent parent
## 47 16
# Mean per group
RT_year <- aggregate(x = f_fledg_mother, by = list(Raising = f_fledg_mother$raising_type, Year = f_fledg_mother$hatch_year),
FUN = mean)## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
GLM
null_model_f_fledg_raising <- glm(formula = nr_f_fledg ~ 1,
data = f_fledg_mother,
family = "poisson")
summary(null_model_f_fledg_raising)##
## Call:
## glm(formula = nr_f_fledg ~ 1, family = "poisson", data = f_fledg_mother)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0541 -1.0541 -1.0541 0.5354 2.2868
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5878 0.1690 -3.478 0.000506 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 70.964 on 62 degrees of freedom
## Residual deviance: 70.964 on 62 degrees of freedom
## AIC: 128.63
##
## Number of Fisher Scoring iterations: 5
f_fledg_raising_glm <- glm(formula = nr_f_fledg ~ raising_type,
data = f_fledg_mother,
family = "poisson")
summary(f_fledg_raising_glm) # non significant##
## Call:
## glm(formula = nr_f_fledg ~ raising_type, family = "poisson",
## data = f_fledg_mother)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1299 -1.1299 -0.7906 0.4177 2.8628
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4490 0.1826 -2.459 0.0139 *
## raising_typeparent -0.7142 0.4830 -1.479 0.1393
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 70.964 on 62 degrees of freedom
## Residual deviance: 68.387 on 61 degrees of freedom
## AIC: 128.05
##
## Number of Fisher Scoring iterations: 6
GLMM
f_fledg_raising_glmm <- glmer(formula = nr_f_fledg ~ raising_type + (1|parent_f),
data = f_fledg_mother,
family = "poisson")
summary(f_fledg_raising_glmm) # not significant## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: poisson ( log )
## Formula: nr_f_fledg ~ raising_type + (1 | parent_f)
## Data: f_fledg_mother
##
## AIC BIC logLik deviance df.resid
## 125.3 131.7 -59.7 119.3 60
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8567 -0.5852 -0.3750 0.3852 2.7081
##
## Random effects:
## Groups Name Variance Std.Dev.
## parent_f (Intercept) 0.6684 0.8176
## Number of obs: 63, groups: parent_f, 31
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8428 0.3509 -2.402 0.0163 *
## raising_typeparent -0.9307 0.7060 -1.318 0.1874
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## rsng_typprn -0.142
Plot
f_fledg_raising_predi <- ggpredict(f_fledg_raising_glmm, terms = "raising_type")
f_fledg_raising_predi## # Predicted counts of nr_f_fledg
## # x = raising_type
##
## x | Predicted | 95% CI
## ----------------------------------------
## foster parent | 0.43 | [0.22, 0.86]
## parent | 0.17 | [0.04, 0.73]
##
## Adjusted for:
## * parent_f = 0 (population-level)
plot(x = f_fledg_raising_predi,
show.title = F,
limits = c(0,1)) +
labs(x = "Raising type",
y = "Number of female fledglings per potential mother",
tags = "") +
theme_classic() +
theme(
axis.line = element_line(colour = "black"),
axis.text = element_text(size=14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
axis.title.x = element_text(colour="black", size = 15),
axis.title.y = element_text(colour="black",size=15, vjust = 5))Testing the models
We use the DHARMa package.
Prepare it
f_fledg_raising_simu <- simulateResiduals(fittedModel = f_fledg_raising_glmm)
plot(f_fledg_raising_simu)Dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.94162, p-value = 0.88
## alternative hypothesis: two.sided
Zero-inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.9154, p-value = 0.488
## alternative hypothesis: two.sided
Temporal Autocorrelation
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 2.1776, p-value = 0.4782
## alternative hypothesis: true autocorrelation is not 0
GLMM of the number of fledglings per potential mother per colony
Data exploration
# Plot the number of fledglings
ggplot(data = f_fledg_mother, mapping = aes(x = nr_f_fledg)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Plot of the number of fledglings per colony per year
f_fledg_mother %>%
group_by(breeding_area, hatch_year) %>%
mutate(sum.fledg = sum(nr_f_fledg)) %>%
ggplot(mapping = aes(x = hatch_year, y = sum.fledg)) +
geom_line(aes(group = breeding_area, color = breeding_area))##
## Burghausen Kuchl allocation
## 28 28 7
Remove individuals with colony “allocation”
f_fledg_mother <- f_fledg_mother[f_fledg_mother$breeding_area != "allocation", ]
f_fledg_mother <- droplevels(f_fledg_mother)
table(f_fledg_mother$breeding_area) # half half##
## Burghausen Kuchl
## 28 28
GLM
null_model_f_fledg_col <- glm(formula = nr_f_fledg ~ 1,
data = f_fledg_mother,
family = "poisson")
summary(null_model_f_fledg_col)##
## Call:
## glm(formula = nr_f_fledg ~ 1, family = "poisson", data = f_fledg_mother)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1180 -1.1180 -1.1180 0.4359 2.1591
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.470 0.169 -2.781 0.00542 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 62.719 on 55 degrees of freedom
## Residual deviance: 62.719 on 55 degrees of freedom
## AIC: 120.39
##
## Number of Fisher Scoring iterations: 5
f_fledg_colony_glm <- glm(formula = nr_f_fledg ~ breeding_area,
data = f_fledg_mother,
family = "poisson")
summary(f_fledg_colony_glm) # not significant##
## Call:
## glm(formula = nr_f_fledg ~ breeding_area, family = "poisson",
## data = f_fledg_mother)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1339 -1.1019 -1.1019 0.4116 2.1909
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.49899 0.24254 -2.057 0.0396 *
## breeding_areaKuchl 0.05716 0.33820 0.169 0.8658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 62.719 on 55 degrees of freedom
## Residual deviance: 62.691 on 54 degrees of freedom
## AIC: 122.36
##
## Number of Fisher Scoring iterations: 6
GLMM
f_fledg_colony_glmm <- glmer(formula = nr_f_fledg ~ breeding_area + (1|parent_f),
data = f_fledg_mother,
family = "poisson")
summary(f_fledg_colony_glmm)## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: poisson ( log )
## Formula: nr_f_fledg ~ breeding_area + (1 | parent_f)
## Data: f_fledg_mother
##
## AIC BIC logLik deviance df.resid
## 121.1 127.2 -57.6 115.1 53
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.9485 -0.5909 -0.5461 0.5472 2.2144
##
## Random effects:
## Groups Name Variance Std.Dev.
## parent_f (Intercept) 0.4809 0.6934
## Number of obs: 56, groups: parent_f, 25
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.77960 0.38967 -2.001 0.0454 *
## breeding_areaKuchl 0.02262 0.47364 0.048 0.9619
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## brdng_rKchl -0.606
Plot
f_fledg_colony_predi <- ggpredict(f_fledg_colony_glmm, terms = "breeding_area")
f_fledg_colony_predi## # Predicted counts of nr_f_fledg
## # x = breeding_area
##
## x | Predicted | 95% CI
## -------------------------------------
## Burghausen | 0.46 | [0.21, 0.98]
## Kuchl | 0.47 | [0.22, 1.01]
##
## Adjusted for:
## * parent_f = 0 (population-level)
plot(x = f_fledg_colony_predi,
show.title = F,
limits = c(0,1.15)) +
labs(x = "Colony",
y = "Number of female fledglings per potential mother",
tags = "") +
theme_classic() +
theme(
axis.line = element_line(colour = "black"),
axis.text = element_text(size=14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
axis.title.x = element_text(colour="black", size = 15),
axis.title.y = element_text(colour="black",size=15, vjust = 5))Testing the models
We use the DHARMa package.
Prepare it
f_fledg_colony_simu <- simulateResiduals(fittedModel = f_fledg_colony_glmm)
plot(f_fledg_colony_simu)Dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.96286, p-value = 0.968
## alternative hypothesis: two.sided
Zero-inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.92394, p-value = 0.608
## alternative hypothesis: two.sided
Temporal Autocorrelation
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.5061, p-value = 0.05923
## alternative hypothesis: true autocorrelation is not 0
Female and male fledglings per nest per year
Based on the values for RRNest ## GLMM of the number of fledglings per nest per raising type
## [1] "b_pair_date" "b_parent_f" "b_parent_f_id" "b_parent_m" "b_parent_m_id" "b_breeding_colony" "b_count_eggs" "b_count_hatching"
## [9] "b_count_fledgling" "breed_year" "raising_type"
Data exploration
# Plot the number of fledglings
ggplot(data = nest_numbers_fm, mapping = aes(x = b_count_fledgling)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Plot of the number of fledglings per raising type per year
nest_numbers_fm %>%
group_by(raising_type, breed_year) %>%
mutate(sum.fledg = sum(b_count_fledgling)) %>%
ggplot(mapping = aes(x = breed_year, y = sum.fledg)) +
geom_line(aes(group = raising_type, color = raising_type))# Sum of the breeding events with reproduction per nest and per raising type
table(nest_numbers_fm$raising_type)##
## foster parent parent added
## 31 3 32
Analysis with temporarily added females as mothers
GLM
null_model_fm_nest_raising_full <- glm(formula = b_count_fledgling ~ 1,
data = nest_numbers_fm,
family = "poisson")
summary(null_model_fm_nest_raising_full)##
## Call:
## glm(formula = b_count_fledgling ~ 1, family = "poisson", data = nest_numbers_fm)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.06706 -0.86863 -0.09432 0.55657 1.66615
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.75911 0.08422 9.014 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 72.857 on 65 degrees of freedom
## Residual deviance: 72.857 on 65 degrees of freedom
## AIC: 227.22
##
## Number of Fisher Scoring iterations: 5
fm_nest_raising_full_glm <- glm(formula = b_count_fledgling ~ raising_type,
data = nest_numbers_fm,
family = "poisson")
summary(fm_nest_raising_full_glm) # non significant##
## Call:
## glm(formula = b_count_fledgling ~ raising_type, family = "poisson",
## data = nest_numbers_fm)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.30940 -0.82416 -0.04514 0.60937 1.27256
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.72490 0.12500 5.799 6.66e-09 ***
## raising_typeparent 0.25593 0.37500 0.682 0.495
## raising_typeadded 0.04347 0.17354 0.251 0.802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 72.857 on 65 degrees of freedom
## Residual deviance: 72.409 on 63 degrees of freedom
## AIC: 230.77
##
## Number of Fisher Scoring iterations: 5
GLMM
fm_nest_raising_full_glmm <- glmer(formula = b_count_fledgling ~ raising_type + (1|b_parent_f),
data = nest_numbers_fm,
family = "poisson")## boundary (singular) fit: see ?isSingular
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: poisson ( log )
## Formula: b_count_fledgling ~ raising_type + (1 | b_parent_f)
## Data: nest_numbers_fm
##
## AIC BIC logLik deviance df.resid
## 232.8 241.5 -112.4 224.8 62
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6330 -0.7409 -0.0449 0.6511 1.4289
##
## Random effects:
## Groups Name Variance Std.Dev.
## b_parent_f (Intercept) 0 0
## Number of obs: 66, groups: b_parent_f, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.72490 0.12500 5.799 6.66e-09 ***
## raising_typeparent 0.25593 0.37500 0.682 0.495
## raising_typeadded 0.04347 0.17354 0.251 0.802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rsng_typp
## rsng_typprn -0.333
## rsng_typddd -0.720 0.240
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Plot
fm_nest_raising_full_predi <- ggpredict(fm_nest_raising_full_glmm, terms = "raising_type")
fm_nest_raising_full_predi## # Predicted counts of b_count_fledgling
## # x = raising_type
##
## x | Predicted | 95% CI
## ----------------------------------------
## foster parent | 2.06 | [1.62, 2.64]
## parent | 2.67 | [1.33, 5.33]
## added | 2.16 | [1.70, 2.73]
##
## Adjusted for:
## * b_parent_f = 0 (population-level)
plot(x = fm_nest_raising_full_predi,
show.title = F,
limits = c(0,6)) +
labs(x = "Raising type",
y = "Number of fledglings per nest (m + f)",
tags = "") +
theme_classic() +
theme(
axis.line = element_line(colour = "black"),
axis.text = element_text(size=14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
axis.title.x = element_text(colour="black", size = 15),
axis.title.y = element_text(colour="black",size=15, vjust = 5))Testing the models
We use the DHARMa package.
Prepare it
fm_nest_raising_full_simu <- simulateResiduals(fittedModel = fm_nest_raising_full_glmm)
plot(fm_nest_raising_full_simu)Dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.89881, p-value = 0.272
## alternative hypothesis: two.sided
Zero-inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.375, p-value = 0.32
## alternative hypothesis: two.sided
Temporal Autocorrelation
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.9305, p-value = 0.7767
## alternative hypothesis: true autocorrelation is not 0
Analysis without temporarily added females as mothers
Here we consider only the raising types foster parent and biological parent
GLM
null_model_fm_nest_raising <- glm(formula = b_count_fledgling ~ 1,
data = nest_numbers_fm_bpfp,
family = "poisson")
summary(null_model_fm_nest_raising)##
## Call:
## glm(formula = b_count_fledgling ~ 1, family = "poisson", data = nest_numbers_fm_bpfp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.05798 -0.66326 -0.08161 0.57021 1.68125
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7503 0.1179 6.367 1.93e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 32.558 on 33 degrees of freedom
## Residual deviance: 32.558 on 33 degrees of freedom
## AIC: 114.34
##
## Number of Fisher Scoring iterations: 5
fm_nest_raising_glm <- glm(formula = b_count_fledgling ~ raising_type,
data = nest_numbers_fm_bpfp,
family = "poisson")
summary(fm_nest_raising_glm) # not significant##
## Call:
## glm(formula = b_count_fledgling ~ raising_type, family = "poisson",
## data = nest_numbers_fm_bpfp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.30940 -0.62940 -0.04514 0.60937 1.27256
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7249 0.1250 5.799 6.66e-09 ***
## raising_typeparent 0.2559 0.3750 0.682 0.495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 32.558 on 33 degrees of freedom
## Residual deviance: 32.122 on 32 degrees of freedom
## AIC: 115.9
##
## Number of Fisher Scoring iterations: 5
GLMM
fm_nest_raising_glmm <- glmer(formula = b_count_fledgling ~ raising_type + (1|b_parent_f),
data = nest_numbers_fm_bpfp,
family = "poisson")## boundary (singular) fit: see ?isSingular
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: poisson ( log )
## Formula: b_count_fledgling ~ raising_type + (1 | b_parent_f)
## Data: nest_numbers_fm_bpfp
##
## AIC BIC logLik deviance df.resid
## 117.9 122.5 -55.9 111.9 31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6330 -0.5669 -0.0449 0.6511 1.4289
##
## Random effects:
## Groups Name Variance Std.Dev.
## b_parent_f (Intercept) 0 0
## Number of obs: 34, groups: b_parent_f, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7249 0.1250 5.799 6.66e-09 ***
## raising_typeparent 0.2559 0.3750 0.682 0.495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## rsng_typprn -0.333
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Plot
fm_nest_raising_predi <- ggpredict(fm_nest_raising_glmm, terms = "raising_type")
fm_nest_raising_predi## # Predicted counts of b_count_fledgling
## # x = raising_type
##
## x | Predicted | 95% CI
## ----------------------------------------
## foster parent | 2.06 | [1.62, 2.64]
## parent | 2.67 | [1.33, 5.33]
##
## Adjusted for:
## * b_parent_f = 0 (population-level)
plot(x = fm_nest_raising_predi,
show.title = F,
limits = c(0,6)) +
labs(x = "Raising type",
y = "Number of fledglings per nest (m + f)",
tags = "") +
theme_classic() +
theme(
axis.line = element_line(colour = "black"),
axis.text = element_text(size=14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
axis.title.x = element_text(colour="black", size = 15),
axis.title.y = element_text(colour="black",size=15, vjust = 5))Testing the models
Testing with the DHARMa package.
Prepare it
fm_nest_raising_simu <- simulateResiduals(fittedModel = fm_nest_raising_glmm)
plot(fm_nest_raising_full_simu)Dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.85891, p-value = 0.312
## alternative hypothesis: two.sided
Zero-inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.2054, p-value = 0.816
## alternative hypothesis: two.sided
Temporal Autocorrelation
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.5941, p-value = 0.2277
## alternative hypothesis: true autocorrelation is not 0
GLMM of the number of fledglings per nest per colony
Data exploration
# Plot the number of fledglings
ggplot(data = nest_numbers_fm, mapping = aes(x = b_count_fledgling)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Plot of the number of fledglings per raising type per year
nest_numbers_fm %>%
group_by(b_breeding_colony, breed_year) %>%
mutate(sum.fledg = sum(b_count_fledgling)) %>%
ggplot(mapping = aes(x = breed_year, y = sum.fledg)) +
geom_line(aes(group = b_breeding_colony, color = b_breeding_colony))# Sum of the breeding events with reproduction per nest and per raising type
table(nest_numbers_fm$b_breeding_colony)##
## Burghausen Kuchl
## 40 26
GLM
null_model_fm_nest_colony <- glm(formula = b_count_fledgling ~ 1,
data = nest_numbers_fm,
family = "poisson")
summary(null_model_fm_nest_colony)##
## Call:
## glm(formula = b_count_fledgling ~ 1, family = "poisson", data = nest_numbers_fm)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.06706 -0.86863 -0.09432 0.55657 1.66615
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.75911 0.08422 9.014 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 72.857 on 65 degrees of freedom
## Residual deviance: 72.857 on 65 degrees of freedom
## AIC: 227.22
##
## Number of Fisher Scoring iterations: 5
fm_nest_colony_glm <- glm(formula = b_count_fledgling ~ b_breeding_colony,
data = nest_numbers_fm,
family = "poisson")
summary(fm_nest_colony_glm) # non significant##
## Call:
## glm(formula = b_count_fledgling ~ b_breeding_colony, family = "poisson",
## data = nest_numbers_fm)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.13037 -0.81506 -0.03506 0.62019 1.73666
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7178 0.1104 6.500 8.02e-11 ***
## b_breeding_colonyKuchl 0.1016 0.1707 0.595 0.552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 72.857 on 65 degrees of freedom
## Residual deviance: 72.504 on 64 degrees of freedom
## AIC: 228.86
##
## Number of Fisher Scoring iterations: 5
GLMM
fm_nest_colony_glmm <- glmer(formula = b_count_fledgling ~ b_breeding_colony + (1|b_parent_f),
data = nest_numbers_fm,
family = "poisson")## boundary (singular) fit: see ?isSingular
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: poisson ( log )
## Formula: b_count_fledgling ~ b_breeding_colony + (1 | b_parent_f)
## Data: nest_numbers_fm
##
## AIC BIC logLik deviance df.resid
## 230.9 237.4 -112.4 224.9 63
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.50640 -0.73335 -0.03492 0.66351 2.06037
##
## Random effects:
## Groups Name Variance Std.Dev.
## b_parent_f (Intercept) 0 0
## Number of obs: 66, groups: b_parent_f, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7178 0.1104 6.500 8.02e-11 ***
## b_breeding_colonyKuchl 0.1016 0.1707 0.595 0.552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## b_brdng_clK -0.647
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Plot
fm_nest_colony_predi <- ggpredict(fm_nest_colony_glmm, terms = "b_breeding_colony")
fm_nest_colony_predi## # Predicted counts of b_count_fledgling
## # x = b_breeding_colony
##
## x | Predicted | 95% CI
## -------------------------------------
## Burghausen | 2.05 | [1.65, 2.55]
## Kuchl | 2.27 | [1.76, 2.93]
##
## Adjusted for:
## * b_parent_f = 0 (population-level)
plot(x = fm_nest_colony_predi,
show.title = F,
limits = c(0,3)) +
labs(x = "Breeding colony",
y = "Number of fledglings per nest (m + f)",
tags = "") +
theme_classic() +
theme(
axis.line = element_line(colour = "black"),
axis.text = element_text(size=14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
axis.title.x = element_text(colour="black", size = 15),
axis.title.y = element_text(colour="black",size=15, vjust = 5))Testing the models
Testing with the DHARMa package.
Prepare it
fm_nest_colony_simu <- simulateResiduals(fittedModel = fm_nest_colony_glmm)
plot(fm_nest_colony_simu)Dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.89981, p-value = 0.304
## alternative hypothesis: two.sided
Zero-inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.3798, p-value = 0.32
## alternative hypothesis: two.sided
Temporal Autocorrelation
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 2.1516, p-value = 0.5357
## alternative hypothesis: true autocorrelation is not 0
Session Info
## DO NOT REMOVE!
## We store the settings of your computer and the current versions of the
## packages used to allow for reproducibility
Sys.time()## [1] "2022-01-06 20:04:48 CET"
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252 LC_MONETARY=German_Germany.1252 LC_NUMERIC=C LC_TIME=German_Germany.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] pwr_1.3-0 lme4_1.1-25 Matrix_1.2-18 ggeffects_1.0.0 effects_4.2-0 carData_3.0-4 DHARMa_0.3.3.0 lubridate_1.7.9 timelineS_0.1.1
## [10] viridisLite_0.3.0 forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2 purrr_0.3.4 readr_1.4.0 tidyr_1.1.2 tibble_3.0.3 tidyverse_1.3.0
## [19] survminer_0.4.8 ggpubr_0.4.0 ggplot2_3.3.2 survival_3.2-7 here_0.1
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_1.4-1 ggsignif_0.6.0 ellipsis_0.3.1 rio_0.5.16 sjlabelled_1.1.7 showtext_0.9 rprojroot_1.3-2 snakecase_0.11.0
## [10] estimability_1.3 fs_1.5.0 rstudioapi_0.11 farver_2.0.3 showtextdb_3.0 remotes_2.2.0 fansi_0.4.1 xml2_1.3.2 codetools_0.2-16
## [19] splines_4.0.3 knitr_1.30 pkgload_1.1.0 jsonlite_1.7.1 nloptr_1.2.2.2 broom_0.7.4 km.ci_0.5-2 dbplyr_1.4.4 compiler_4.0.3
## [28] d6_0.1.0.0 httr_1.4.2 backports_1.1.10 assertthat_0.2.1 survey_4.0 cli_2.0.2 htmltools_0.5.0 prettyunits_1.1.1 tools_4.0.3
## [37] gtable_0.3.0 glue_1.4.2 tinytex_0.26 Rcpp_1.0.5 cellranger_1.1.0 vctrs_0.3.4 nlme_3.1-149 iterators_1.0.13 lmtest_0.9-38
## [46] insight_0.11.1 xfun_0.18 ps_1.4.0 openxlsx_4.2.3 testthat_2.3.2 rvest_0.3.6 lifecycle_0.2.0 devtools_2.3.2 statmod_1.4.35
## [55] rstatix_0.6.0 MASS_7.3-53 zoo_1.8-8 scales_1.1.1 hms_0.5.3 yaml_2.2.1 curl_4.3 memoise_1.1.0 gridExtra_2.3
## [64] KMsurv_0.1-5 stringi_1.5.3 desc_1.2.0 gap_1.2.2 foreach_1.5.1 boot_1.3-25 pkgbuild_1.1.0 zip_2.1.1 rlang_0.4.7
## [73] pkgconfig_2.0.3 evaluate_0.14 lattice_0.20-41 labeling_0.3 processx_3.4.4 tidyselect_1.1.0 magrittr_1.5 bookdown_0.20 R6_2.4.1
## [82] generics_0.0.2 DBI_1.1.0 pillar_1.4.6 haven_2.3.1 foreign_0.8-80 withr_2.3.0 nnet_7.3-14 abind_1.4-5 modelr_0.1.8
## [91] crayon_1.3.4 car_3.0-10 survMisc_0.5.5 rmarkdown_2.4 sysfonts_0.8.1 usethis_1.6.3 grid_4.0.3 readxl_1.3.1 data.table_1.13.2
## [100] blob_1.2.1 callr_3.5.0 rmdformats_0.3.7 reprex_0.3.0 digest_0.6.25 xtable_1.8-4 munsell_0.5.0 mitools_2.4 sessioninfo_1.1.1